##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone/GCClone_gene.all_record.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/DriverChoose/GeneLOH"
    type <- "IGC"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file
type <- opt$type

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info <- data.frame(fread(sample_info))
dat_gene <- data.frame(fread(gene_list))
dat_ccf <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
## 关注的病理亚型
dat_info <- subset( dat_info , Type != "IM + IGC + DGC")
dat_ccf <- subset( dat_ccf , ID %in% dat_info$ID )

dat_ccf <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol %in% dat_gene$Gene_Symbol )

###########################################################################################
## 共享的突变占比
dat_ccf_share <- subset( dat_ccf , Class != "IM" & Share == TRUE )
dat_ccf_share <- data.frame(table(dat_ccf_share[,c("Hugo_Symbol" , "Class")]))
dat_ccf_share$Type <- "Share"

## 非共享的突变占比
dat_ccf_noshare <- subset( dat_ccf , Class != "IM" & Share == FALSE )
dat_ccf_noshare <- data.frame(table(dat_ccf_noshare[,c("Hugo_Symbol" , "Class")]))
dat_ccf_noshare$Type <- "Private"

## 合并
result_plot <- rbind( dat_ccf_noshare , dat_ccf_share )

## 计算总的数量
result_plot <- result_plot %>%
group_by( Hugo_Symbol , Class ) %>%
summarize( Type = Type , Freq = Freq , total_num = sum(Freq) )

result_plot$Ratio <- result_plot$Freq/result_plot$total_num
result_plot$Ratio[is.na(result_plot$Ratio)] <- 0
result_plot$value_text <- paste0( round(result_plot$Ratio , 2) * 100 , "%")

###########################################################################################
## 计算p值

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

if(1!=1){
    dat <- result_plot
    result <- c()
    for( geneN in unique(dat$Hugo_Symbol) ){
        tmp <- subset( dat , Hugo_Symbol == geneN )

        a <- subset( tmp , Class == "DGC" & Type == "Private" )$Freq
        b <- subset( tmp , Class == "DGC" & Type == "Share" )$Freq
        c <- subset( tmp , Class == "IGC" & Type == "Private" )$Freq
        d <- subset( tmp , Class == "IGC" & Type == "Share" )$Freq

        if(length(a)==0){a=0}
        if(length(b)==0){b=0}
        if(length(c)==0){c=0}
        if(length(d)==0){d=0}

        p <- fisher.test( matrix( c(a,b,c,d) , ncol = 2 ) )$p.value

        if( p < 0.001 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
        }

        tmp$p_text <- ""
        tmp$p_text[1] <- p_text
        result <- rbind( result , tmp )
    } 
}


show_gene <- c("TP53" , "APC" , "PIK3CA" , "CDH1")
result <- subset( result_plot , Hugo_Symbol %in% show_gene  )
#gene_order <- c(
#    "TP53" , "ARID1A" , "CDH1" , "APC" , 
#    "ERBB2" , "PIK3CA" , "RNF43" , "MAP2K7" ,
#    "MTRR" , "MUC6" , "CFTR" , "BMP6" , "GAL3ST3")

result$Hugo_Symbol <- factor(result$Hugo_Symbol , levels = show_gene )

plot <- ggplot( data = result , aes( x = Var1 , y = Ratio , fill = Var2 ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    facet_grid(.~Hugo_Symbol) +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".LOH.pdf")
ggsave( out_name , plot , width = 3 , height = 4 )
